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We investigate the dynamics of the critical isoscalar condensate, formed during heavy-ion 
collisions. Our analysis is based on a simplified model where the sigma and the pions are 
the only degrees of freedom. In field description, both in physical and momentum space, 
we find that the freeze-out profile presents a structure which reveals clear traces of the 
critical fluctuations in the sigma-component. In particle representation, using Monte-Carlo 
simulations and factorial moment analysis, we show that signatures of the initial criticality 
survive at the detected pions. We propose the distribution of suitably defined intermittency 
indices, incorporating dynamical effects due to sigma-pion interaction, as the basic observable 
for the exploration of critical fluctuations in heavy-ion collision experiments. 

I. INTRODUCTION 

Experiments of a new generation, with relativistic nuclei, at RHIC and SPS are currently under 
consideration with the aim to intensify the search for the existence and location of the QCD critical 
point in the phase diagram of strongly interacting matter [1, 2]. Important developments in lattice 
QCD [3] and studies of hadronic matter at high temperatures [4] suggest that the QCD critical 
endpoint is likely to be located within reach at SPS energies. A decisive observation, in these 
experiments, associated with the development of a second-order phase transition at the critical 
endpoint of QCD matter is the establishment of power laws in momentum space (self-similarity) 
in close analogy to the phenomenon of critical opalescence in QED matter [5] . These power laws 
reflect the fractal geometry in real space, and a characteristic index of the critical behavior is the 
fractal mass dimension df2 which measures the strength of the order-parameter fluctuations within 
the universality class of critical QCD [6, 7]. 

The physics underlying the endpoint singularity in the QCD phase diagram is associated with 
the phenomenon of chiral phase transition, a fundamental property of strong interactions in the 
limit of zero quark masses. In this case and for a given chemical potential fj,B there exists a critical 
temperature T cr above which chiral symmetry is restored and as temperature decreases below T cr , 
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the system enters into the chirally broken phase of observable hadrons. It is believed that, for two 
flavors and zero quark masses, there is a first-order phase transition line on the (T,hb) surface 
at large \xb [8]. This line ends at a tricritical point beyond which the phase transitions become 
of second order. In the case of real QCD with non-zero quark masses, chiral symmetry is broken 
explicitly and the first-order line ends at a critical point beyond which the second-order transitions 
are replaced by analytical crossovers [9]. 

The order parameter of this transition is the chiral field $ = (a, vr) formed by the scalar, 
isoscalar field a together with the pseudoscalar, isovector field tt = (tt + , 7r°, 7T~). Both fields are 
massless at the tricritical point and when the symmetry is restored at high temperatures, their 
expectation value vanishes, (tt) = (a) = 0. However, in the presence of an explicit symmetry 
breaking mechanism (non-zero quark masses) the pion and sigma fields are disentangled at the 
level of the order parameter of the QCD critical point, which is now formed by the sigma field 
alone. The expectation value of the u-field remains small but not zero near the critical point so 
that the chiral symmetry is never completely restored. The valuable observables in this case are 
associated with the fluctuations of the c-field, (5a) 2 ~ (c 2 ), and they incorporate, in principle, the 
singular behavior of baryon-number susceptibility and in particular the power-law behavior of the 
(T-field correlator [10]. 

The aim of this work is to study the evolution of critical correlations during the development 
of a collision of heavy ions (A + A), in the neighborhood of the QCD critical point. In the initial 
state we assume that the system has reached the critical point in thermal equilibrium so that 
the cr-field fluctuations are described by the 3D Ising critical exponent 5 and in particular by the 
fractal dimension dj2 = (5 ~ 5). The crucial question from the observational point of view 

is whether traces of the initial criticality can survive in the freeze-out regime, which follows the 
equilibration stage, and leave their imprints at the detectors. In [11] we showed that this is indeed 
the case for a small class of events. However the treatment presented in [11] has two important 
disadvantages: (i) the analysis of the fluctuations has been performed in the experimentally inac- 
cessible configuration space and (ii) the determination of the small class of events carrying critical 
characteristics is by itself a non trivial task. In the present investigation we substantially improve 
our approach in this direction, achieving at the same time two main goals: (i) better adaptation 
to experimental requirements, by considering the evolution of the geometry of the critical conden- 
sate through intermittency analysis in transverse momentum space and (ii) definition of suitable 
observables which allow for a quantitative description of the deformation of the critical system by 
considering all the events in the ensemble. To proceed in this direction we introduce the concept 
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of evolutionary intermittency in transverse momentum space, which is detectable in heavy-ion col- 
lision experiments through factorial moment analysis in an event-by-event basis. Independently 
from the duration of the freeze-out process, the final profile of the deformed system presents a 
structure that incorporates all its preceding dynamical history. For its analysis we investigate the 
distribution of the intermittency indices calculated during the evolution towards the final state. It 
is shown that this distribution is sensitive to the initial critical state offering a, free from statistical 
restrictions, tool to explore critical fluctuations. 

The dynamics of the system is fixed by a two-field Lagrangian, C(a, tt), together with appropri- 
ate initial conditions. The out-of-equilibrium phenomena are generated by the exchange of energy 
between the <7-field and the environment which consists of massive pions initially in thermal equi- 
librium. We adopt the picture of a rapid expansion (quench) which is a realistic possibility in this 
experimental framework. We study the out-of-equilibrium evolution of the cr-field and we impose 
its decay to pions according to the characteristics of the a — ► 7r + 7r~ procedure. Since in every 
collisional event the u-particles decay independently, time integrations must be taken into account. 
The freeze-out profile obtained in this way is analyzed both in the field description and particle 
representation. In order to acquire the particle picture, which is the one obtained in experiments, 
we use Monte-Carlo simulations of the aforementioned procedure, and the fractal measures are 
quantified through factorial moment analysis [12, 13]. 

In section II the formulation of the problem and in particular the equations of motion and the 
initial conditions are presented. In section III numerical solutions are given in the context of field 
description in physical and momentum space whereas in section IV we present numerical results 
in particle framework. Finally, in V our results and conclusions are summarized. 



In our approach we assume an initial critical state of the system in thermal equilibrium, dis- 
turbed by a two-field potential V(cr, tt), in an effective description inspired by the chiral theory of 
strong interactions [14, 15]. The Lagrangian density is 



II. THE MODEL 




(1) 



with the potential 
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where a = a(x,t) and 7? = w(x,t). The potential has the usual a- model form plus a term which 
breaks the symmetry along the cr-direction. With the addition of the mass term for the pion field, 
we ensure that it has a constant mass equal to m n . Finally, the constant terms in (2) shift the 
value of the potential at the minimum to zero. We fix the parameters of the Lagrangian using the 
phenomenological values m,r~ 139 MeV, vq ~ 87.4 MeV and A 2 ~ 20 [16]. 
The equations of motion resulting from (1) are: 

a - V 2 cr + X 2 (a 2 + tt 2 - vl)a = v ?nl 

V 2 w + A 2 (o- 2 + 7f 2 -^)vf + m 2 vr = 0, (3) 

where vr 2 = (vr+) 2 + (tt ) 2 + (-ir~) 2 . 

Using a constant value vq in eq. (2) implies a non-vanishing mass (finite correlation length) for 



the cr-field, m a = ^J2\ 2 Vq, in contradiction with its critical profile. This inconsistency is restored if 
we assume a finite-time mechanism instead of the instant quench i.e. the instantaneous formation 
of the potential (2) at t = 0. The simplest model which leaves the equations of motion unaffected, 
used also in cosmological phase transitions [17], is the so-called linear quench [11, IS]. It assumes 
that the minimum v of the potential increases linearly with time, starting from zero and ending at 
the zero temperature value vq ~ 87.4 MeV after a time interval r which is the quench duration: 

{vnt/r for t < t 

(4) 
vq for t > r. 

This assumption leads naturally to a known evolution of m a (t) = d 2 V(a, i:) / da 2 \ a=(Jmin ^ = j? rnin : 

{ ^ for t < t 

m*(t) = { ; T ~ (5) 

[ V2>?vo for t > t. 

Thus, m a presents the expected behavior, i.e. it is m a = at t = (critical cr-field) and, with 
increasing time, approaches and acquires its zero temperature value m a ^ = \^2vqX. 

We investigate the evolution of the above system using initial field configurations dictated by 
the onset of the critical behavior. In this case we expect that the cr-field, being the order parameter, 
will possess critical fluctuations, whereas the 7r-fields are thermal, and the entire system will be 
in thermodynamical and chemical equilibrium. Obviously, the subsequent evolution, determined 
by eqs. (3), will generate strong deviations from equilibrium. The question is whether the initial 
critical profile (presenting the critical point information) survives for times large enough in order 
to leave imprints at the freezeout. Before going on with the detailed study of the dynamics, we 
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first refer to the generation of an ensemble of critical field configurations (initial conditions), and 
we provide the necessary tools for its quantitative description. 

The detailed procedure for the production of the critical ensemble is given in [19]. Interpreting 
the square of the a-field as local density, we use the effective action at the critical point [20]: 

T[a}= [ d D xd(Va) 2 +ga 5+1 }, (6) 

JR 1 

(where 5 is the isothermal critical exponent of the corresponding universality class and g is an 
associated effective coupling) in order to weight the sum over field configurations in the partition 
function. This sum is saturated through saddle point solutions which consist approximately of 
piecewise constant configurations extended over domains of variable size [21]. In this way we 
record a large number of statistically independent configurations possessing critical characteristics, 
and the corresponding behavior is described by a fractal measure demonstrated in the power-law 
dependence of the mean "mass" m 2x (R) on the distance R around a point xq defined by: 

m 2x (R) = [ |cj 2 (f)a 2 (f )| d D xd D x ~ R d f 2 *, (7) 

JR 

for every xq. dt 2x is the fractal mass dimension of a specific configuration [22, 23], and the subscript 
2 marks the use of a 2 for its definition, while the subscript x indicates the calculation in physical 
space. df2x is related to space dimensionality D and isothermal critical exponent 5 through: 

_ D(5-l) 

d ^ ~ -8 + T- (8) 

For the 3D Ising universality class, g ~ 2, D = 3, 5 ~ 5, therefore df2 X \D=3 ~ 2. The 3D 
Ising description has to be adapted to the physical framework of colliding nuclei. The considered 
system may be affected by non-equilibrium processes along the collision axis, and therefore the 
picture of statistical equilibration is valid safely only for the transverse 2D section. Moreover, the 
relativistic nature of the longitudinal motion needs an extra assumption (Feynman- Wilson fluid 
hypothesis [24]) in order to formulate global thermal equilibrium in terms of space-time rapidity. 
It is then natural to consider the 2D projection of the original 3D Ising system as a description 
of the fireball in the transverse space. Thus, assuming a cylindrically symmetric fractal, its 2D 
projection has a mass dimension df 2x \2D proj = df 2x ~ 4/3, and the a-field depends only on the 
transverse coordinates. We produce an ensemble of 10 4 cr-configurations, each one possessing its 
own fractal mass dimension, and these df 2x ^s satisfy (7) and (8) within an error of 2%. 
Equivalently, we can define the fractal mass dimension in (transverse) fc-space as: 

m 2k (K)= [ \d 2 (k)d 2 {k Q )\d 2 kd 2 k ~K d v\ (9) 

JK 



where ^(fc) is the 2D Fourier transformation of cr 2 (x), and K is the "distance" around a point ko 
in the transverse momentum space. Obviously, df2 X and df2k are related through: 

df2k = 2-df 2x - (10) 

Performing numerically the aforementioned Fourier transformation, we result in a critical ensemble 
in A;-space satisfying (9), whereas (10) holds within an error of less than 10 -2 per cent. 

Having established the production algorithm for the critical ensemble of the cr-field, and the 
necessary measures for its quantitative description in physical and fc-space, we can proceed to its 
time evolution. As we have already mentioned in the introduction, we are interested in inves- 
tigating the evolution of the system determined by equations (3), using as initial conditions an 
ensemble of 10 4 independent cr-configurations on the lattice possessing critical characteristics, and 
10 4 configurations, for each ir component, corresponding to an ideal gas at temperature To ~ 140 
MeV, which simulates the pion environment in a heavy-ion collision. The production of thermal 
7r-configurations has been described in [11, 25]. 

III. NUMERICAL SOLUTIONS: PHYSICAL AND /c-SPACE DESCRIPTION 

We evolve the critical system according to (3) and we determine the field values at any time. The 
quantitative deformation of the initial fractal geometry, is embedded in the evolution of m,2x(R) 
(or of m,2k(K) in fc-space) and in particular in the evolution of the exponent of its power-law 
dependence on R. We call this exponent ip x {t) (similarly ipk{t) m A;-space). Initially ip x (0) = df2 X 
(and ipkifi) = df2k), i-e. it is the fractal mass dimension of the critical system. 

As we showed in [11] ijj x (t) reaches in short time interval the value of the embedding dimension, 
which in the present case, due to the projection into transverse space, is equal to 2, and subsequently 
fluctuates becoming almost equal to ip x (0) at particular times. Thus, the initial critical profile 
reappears partially and quite periodically at times when the spatial mean value of the cr-field 
returns close to its initial value (which is (cr(0)) ~ as implied by the critical behavior). Finally, 
after some oscillations ip x {t) relaxes very close to 2, i.e. the initial fractal geometry is completely 
lost. The detailed explanation of this phenomenon is given in [26]. 

In fig. 1 we depict the evolution of ip x (t) and ipkit) f° r ° ne configuration of the ensemble, using 
for the quench duration the value r = 5 fm, where the calculations are achieved through the 
integrals (7) and (9) at every time. We observe the partial re-establishment of the initial critical 
profile, and its final complete deformation. Furthermore, we verify that relation (10) holds in very 
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good precision at all times. Lastly, a comment must be made for the overflow of ip x {t) above the 
embedding dimension value 2, and equivalently for the negatives values of ipk(t)- Actually, this 
result is artificial and it arises from small disturbance, due to finite-size effects, of the linearity 
of In m2a; (R) versus In R, at these time intervals. However, if we proceed to larger lattices, the 
exceeding of 2 becomes less strong, as we have shown in [26]. This is possible for one event, but 
it demands extremely huge computational times for the whole ensemble. Therefore, we remain in 
the aforementioned selection, having in mind that slopes larger than 2 (or less than 0) should be 
consider as being equal to 2 (or 0). Alternatively, we can avoid this overflow by enforcing a cut 
on the fit accuracy. As in this work we are interested in the reappearance of the initial fractal 
geometry, our results are trustworthy being not affected by this behavior. 

Now, let us examine the evolution of the whole critical ensemble. To keep the information 
of the evolution of each individual configuration, i.e. we perform an event-by-event analysis of 
the ensemble evolution. The reason is that, as we showed in [11], the quantitative characteristics 
of the initial critical behavior sustaining are more transparent in fluctuation-related measures 
rather than in ensemble-averaged ones. Thus, at each time instant we calculate a histogram for 
the distribution p(ip x ), which corresponds to the individual values of ifi x (t) for each of the 10 4 
independent ^-configurations. In fig. 2 we depict p{ip x ) for six successive times. As we can see, 
initially all the <r-configurations have ip x (ty — df2x equal to the theoretical value 4/3 (arising from 
(8) with D = 2, 5 = 5) within an error of 2%. As time passes, ip x {t) for each configuration evolves 
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FIG. 2: (color online) Distribution p(ip x ) for the individual values of ^p x (t) for each of the 10 4 independent 
a -configurations, for six successive times, for quench duration r = 5 fm. 



independently (similarly to fig. 1) and both the maximum and the standard deviation of p(ip x ) 
changes. However, as we can clearly observe, the non-trivial fractal geometry prevails for some 
time in a subset of the ensemble, and eventually disappears leading to p{i^ x ) centered at the value of 
the embedding dimension 2, with decreasing standard deviation, i.e. all configurations do gradually 
lose their initial critical profile. Similarly, we repeat the same steps for the /c-space-defined ipk(t), 
and the corresponding evolution of p(ipk) is depicted in fig. 3. The explanation of the behavior of 
p(ipk) is straightforward, and as expected the loss of the initial critical behavior is reflected to an 
increasingly narrower distribution p(ipk) centered at zero. 

The above evolution of independent configurations is unrealistic and is only useful in order to 
demonstrate the fractal geometry reappearance. In fact, after m a (t) reaches the threshold of 2m n 
at time tth which according to (5) is: 

the ex-particles can decay to pions through a — ► ir + Tr~ . Thus, traces of the initial criticality can be 
transferred to the produced pions. In [11] we were interested mainly in the qualitative description 
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FIG. 3: (color online) Distribution p(ipk) for the individual values of ipk(t) for each of the 10 4 independent 
a -configurations, for six successive times, for quench duration r = 5 fm. 



of this phenomenon and we showed that signs of the initial critical behavior can indeed be found 
in the pion sector. In the present work we proceed to a detailed and more quantitative analysis. 
To achieve this, we first have to specify the width of the cr-decay. This is given by the expression: 



r(t) = Toy 1 — (2m 7r /m cr (i)) 2 , which is easily treated in our approach since we have assumed a 
specific (and quite general) form for the quench evolution and consequently for m a (t) (relation 
(5)). Therefore, for the decay width we acquire: 



r(t) 



for t < t th 



Wi-^OT for hh<t<r (12) 



V s 



W 1 -^ for T<t - 



The factor To can be determined by fixing the value of the c-width at zero temperature. Eq. (12) 
allows us to determine the probability density pdec{t) that a a-particle will decay at t, using the 
definition of T(t) in the particle picture: ^jf- = —T(t)N a and neglecting any influence of the 
surrounding thermal medium, as well as regeneration effects. The probability Pit) to decay a 
sigma particle in the time interval [0, t), assuming that each sigma decays statistically independent 
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from the others, is then given as: 



leading to: 



jy,(o)-jy g (t) rm dz 

{) N a (0) 



We finally get the expression: 



Pdec(t) 







for t < t th 



Amir 2 



-To 



< 2_±^ T 2_2 Ir ^ Tarccos( 2m 2L r ) 



(7,0 



ct.O* 



4 "4 p 



-Tot J I 



^2^r rarccos(^2^) 
m a,0 m cr,0 



for 



for 



tfh <t<T 



T <t. 



(13) 



Equation (13) leads to a distribution for the mass of the decaying cr's with a characteristic peak 
at ml, which for r — ► oo tends to the threshold value 2m n ( tti^j | r — — ^ 27?!^). Qualitatively, 



this behavior corresponds to the singular enhancement of the spectral density of the sigma field 



in medium (p a {u) 



1 



but 



) in the case of partial symmetry restoration [27]. In practice, 
using r = 5 fm we get m* ~ 282 MeV (leading to a decay width ~ 130 MeV) which is very close 
to the two-pion threshold, and the associated cr-mass distribution resembles, within our simplified 
approach, the in- medium spectral properties of the sigma field at the critical point. In the upper 
graph of fig. 4 we depict T(t) for r = 5 fm, taking To = 800 MeV and m CT) o ~ 560 MeV, while in 
the lower graph we show the corresponding Pdec{t)- 

Having formulated the quantitative treatment of a-decay we can proceed to our investigation. 
The scenario we explore is the following: We evolve the ensemble of 10 4 critical ^-configurations in 
an environment of thermal pions. After tth each one of them will decay independently according 
to Pdec{t) given by (13), i.e. each configuration will decay at each own t^eo an d this corresponds 
to a single event. However, experimentally we get only a freeze-out profile of events, since there is 
no available information concerning the decay time at which the pions, observed at the detectors, 
were produced. To acquire this profile we have to dress the evolution of pfyx) with the decay 
probability Pdecif)- Therefore, the event-by-event freeze-out profile in physical space will be the 
integral of p^ x (t) weighted with pdec(t), up to a final time corresponding to the freeze-out time 

tfrz • 

ftfrz 

- / 
JO 



Pwip x (tfr 



P^ x (t)Pdec(t)dt. 



(14) 
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FIG. 4: Decay width T(t) and p d ec(t) for r = 5 fm, T = 800 MeV and m a . w 560 MeV. 



Although in general the weighted integral has to be performed up to infinite freeze-out time, 
(or up to its experimentally determined value), due to the exponential decrease of Pdecif) all the 
profiles with freeze-out times tf rz larger than t* ~ 10/m are practically identical. Similarly, the 
corresponding profile in fc-space will be the weighted integral p w ^ k (tf rz ) of p^ k {t) with Pdec(t)- 

In fig. 5 we depict this freeze-out profile, both in physical and /c-space, for freeze-out time 
tf rz = 10 fm (normalized to give the correct event-multiplicity, i.e. 10 4 ). As we observe, apart 
from the main peak at 2 in physical space (the main peak at in fc-space) which corresponds to 
complete loss of the initial critical behavior, there is a secondary "shoulder" which is a result of the 
time history of all the events that preserved traces of their initial critical profile at their decay-time 
tdec- This secondary structure is a safe and robust signature of the initial criticality. Furthermore, 
it is also univocal, since starting from conventional (non-critical) initial conditions it is impossible 
to acquire even one non-conventional event at freeze-out, as we have shown in [11]. 

IV. NUMERICAL SOLUTIONS: MONTE-CARLO SIMULATION 

In the previous section we investigated the evolution of the critical system, constituted by a 
large number of individual events, based on the elaboration of the a and n fields. In this section 
we proceed to a more advanced description and refer to particles, since these are the entities that 
are detectable. This process demands a sequence of steps in order to bring our treatment closer to 
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FIG. 5: On the left the freeze-out profile in physical space, given as the weighted integral Pwip^itfrz) of 
Pip* W with Pdec{t), for t f rz = 10 fm. On the right the corresponding profile in k-space, given as the weighted 
integral p w ^ k (tf rz ) of p^ k (t) with Pdec(f), for the same freeze-out time. 



the conditions present in a heavy- ion collision experiment. 

In our study we use cx-configurations which correspond to single events and we assign to each 
of them a ^-multiplicity of ~ 140 er's. This corresponds to A + A collisions with intermediate 
size nuclei assuming the ideal case when the charged pions observed in the (hypothetical) detector 
originate exclusively from the decay of these sigmas and they do not suffer final state interactions. In 
fact, the integral: f v a 2 (x)dx is much larger than the integral: f v 7r 2 (x)dx for all times. Therefore, 
interpreting a 2 as the sigma-particle density and 7r 2 as the pion density, the above property indicates 
that at criticality the multiplicity of the pions in the thermal environment is much smaller than 
the corresponding cr-multiplicity. Now, the presence of a large number of sigma-particles per event 
can dramatically complicate things since every c-particle within an event can decay independently, 
and there is no consistent formalism for treating the subsequent evolution of the u-field after its 
partial decay. However, since Pdec(t) has a sharp peak, we assume that all the cr's within a single 
event decay at the same tdec, which is of course varying for the different events of the ensemble. 

Conclusively we follow this procedure: We evolve each configuration (individual event) according 
to (3). We determine the corresponding tdec according to Pdec(t), and at this specific time we find 
the 2D Fourier transformation &2(k) of the square of the cr-field a 2 (x). Then we produce the 
momentum components of each of the 140 particles through a Monte-Carlo simulation of this 
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momentum profile. Knowing the momenta of the cr-particles we can easily use the kinematics of 
the decay a — > tt + it~ , in order to determine the momenta of the produced pions, similarly to the 
previous work [7]. Note that our algorithm allows a loss of information through the secondary 
decay a — ► 2ir°, and therefore the produced pion pairs are less than the 140 <r's (they vary between 
80 and 100, i.e. around the two thirds of the cr-multiplicity) . Finally, we repeat this procedure for 
the whole ensemble of events. 

In particle representation, the requisite statistical tool for the quantitative description of the 
fractal properties is the factorial moment analysis. We first perform a partition of the two- 
dimensional transverse fc-space into M 2 cells, and we c&lculcite the particle multiplicity inside 
each cell. Then we repeat this for all possible cell sizes (between two scales). The corresponding 
second factorial moment is defined as: 



and for a critical profile, corresponding to isothermal critical exponent 5, this two-dimensional 
transverse moment analysis leads to the power law: 



with the exponent /2 given by: 



ji is called "intermittency index", since (16) determines the self-similarity of the multi-particle 
state spectrum, when it is being analyzed in different scales, the so-called intermittency [7, 12, 13]. 
Note that since we work in the framework of event-by-event analysis throughout this study, we 
define the second factorial moment for each of the individual configurations, and thus we do not 
perform event averaging in (15) which is usually taken in the literature. This leads to a slight 
decrease in the quality of the fit in i^(M) vs M 2 graph, due to reduced statistics. However, as we 
have already mentioned, the event-by-event analysis reveals rich information about the fluctuations 
evolution, which remains hidden in an average analysis. 

There is one last problem that remains to be solved, before the framework of our approach is 
completed. Although, given the a-particle momenta, we can verify relation (16) in large accuracy, 
the pion analysis presents a complication. In an event of a heavy-ion collision the majority of 
the cr-particles decays to a pair of ir + ,ir~ , while the rest decay to neutral pions which cannot 
be detected straightforward. However, in the detectors we record all the charged pions of an 




(15) 




(16) 
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event simultaneously, and thus we cannot find which pairs of them come from the same a, in 
order to reconstruct them and compute the factorial moments correctly. A more sophisticated 
reconstruction procedure is necessary, and it has been elaborated in [6]. The basic features are 
the following: We first choose a narrow invariant mass window, located just above the two-pion 
threshold, to perform the reconstruction of the sigma field. This is dictated by the form of the 
mo-distribution (which for r = 5 fm is peaked at m* ~ 282 MeV) and is compatible with the 
spectral properties of the critical sigma-field in medium due to the partial restoration of chiral 
symmetry [27]. We then combine all n + with all ir~ within a single event making di-pions, and 
we record the subset of them that has invariant mass inside the selected window. Obviously, 
some of these di-pions correspond to real cr's, but some of them (usually the majority) consist a 
combinatorial background modifying the value of the factorial moments. This background can be 
simulated by di-pions reconstructed from mixed events and it is removed by the substraction of the 
corresponding moments. The corrected factorial moments reveal, to a large extent, the behavior of 
the cr's before their decay. Moreover, this substraction eliminates also the thermal pions that are 
present in our model and give conventional contribution to the factorial moments. We point here 
that in the reconstruction procedure there is an additional information loss, concerning the initial 
sigma state, in the neutral pion channel. Furthermore, performing it event by event it stiffens the 
procedure due to the low statistics. However, in the end we re-acquire an exponent of the factorial 
moments that is equal to the expected within an error of less than 5%. 

Having formulated our method, we can proceed to numerical simulation, verifying first the va- 
lidity of our approach at each step. In order to test the Monte-Carlo procedure and the factorial 
moment analysis in the a-sector, we simulate the decay of the initial critical <r-configurations (cor- 
responding to 8 = 5) and we perform moment analysis on the produced cr-particles. In subfigure 6a 
we depict the second factorial moment for each one of the events and in subfigure 6b we present 
the corresponding R 2 values which control the quality of the linear fit. R 2 is defined as: 

r2 = (<yy> - <yxy >) 2 

(< y 2 > - < y > 2 )(< y 2 > - < y > 2 ) ' 

where y are the values for a given measure Y while y are the corresponding values of the fitting 
function. Each value of these histograms arise from a single event, where the corresponding second 
factorial moment plot is like the one shown in subfigure 6c (here we show an event with slope 
that is one fifth standard deviations from the histogram average). As we can see from the p{$2) 
histogram, the obtained fi values in the a-particle sector are distributed around the theoretical 
expectation /2 = 2/3 according to relation (17). The width of the distribution in general depends 
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FIG. 6: (color online)/^ graph a) we present the histogram p{$2) of the exponents of the second factorial 
moment for each one of the events, in the case of a -particles produced directly from the critical ensemble 
corresponding to 5 = 5. In graph b) we show the corresponding values of R 2 which determine the precision 
of the fits. In graph c) we depict the second factorial moment of a single event ( with slope less than one fifth 
standard deviations away from the histogram average). In graph d) we present the average second factorial 
moment, which is usually used in the literature. 



on the multiplicity of c-particles per event which in our case is fixed to 140. Thus, within our 
treatment the critical profile of the cr-particles, through both Monte-Carlo simulation and factorial 
moment analysis, is reproduced in a consistent way. Lastly, in fig. 6d we depict the usually 
calculated average second factorial moment. Although linearity (described by R 2 ) is here better 
and the obtained slope is close to the expected one, the average calculation has the disadvantage of 
information loss due to the replacement of the entire structure of the p(/2)-histogram by a single 
value. As we see in the following the overall structure of the intermittency index histogram proves 
to be very important and this is the reason why we insist on the event-by-event analysis. Finally, 
note here that the average factorial moment of fig. 6d is not exactly the average of the histogram 
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p{f 2)1 since in the former case the mean is calculated before the sum in relation (15), while in the 
latter we take the average of the slopes. 

Having verified our approach at the level of the initial a-configuration we can proceed to the time 
evolution. We evolve the system of 10 4 critical a-configurations (events) in an environment of 10 4 
thermal pion-configurations and we let them decay independently according to Pdec(t). Prom each 
configuration we get 140 cr-particle momenta through Monte-Carlo simulation, and using kinematics 
we perform the decay with a branching ratio 2:3 to 7r + ,7r~ pairs. These are simultaneously detected 
at the freeze-out (similarly to the previous section we use a value of tt rz = 10 fm which is by far 
satisfactory), prohibiting the identification of the original cr's. Thus, as stated above, due to the 
combinatorial background, a simple factorial moment analysis will not be sufficient to reveal the 
critical properties of the decaying system. Indeed, in fig. 7 we present the corresponding four 
graphs of the moment analysis at the freeze-out. In subfigure 7a we show the histogram p(Q of 
the exponent £ for the charged pions, which for a fractal profile is equal to /2 = 2/3 according to 
(17), and in subfig. 7b we give the corresponding R 2 values. In subfig. 7c we depict the second 
factorial moment for a single event (with slope less than one third standard deviations away from 
the histogram average) while in 7d we present the average factorial moment for the whole ensemble 
of events. As we mentioned, the simultaneous detection of all the pions suppresses in this sector the 
signal associated with the partial restoration of the critical sigma profile, and reduces the obtained 
£'s values. Thus, in order to recover the signature of the initial fractal geometry, the previously 
described more sophisticated reconstruction procedure is necessary. 

In fig. 8 we depict the results of the factorial moment analysis for the pions at the freeze-out, 
after the reconstruction step, using an analysis window of 10 MeV just above the two-pion thresh- 
old 2m n . This particular choice is appropriate since it includes the ra^+^~ invariant mass value 
which corresponds to the maximum of the mass distribution of the decaying sigmas (m* ~ 282 
MeV). In addition, in this kinematical window the probabilities of the appearance of (-7r + 7r~)-pairs 
originating from the decay of conventional resonances (as for example the p-meson) is vanishingly 
small. Furthermore, performing our analysis so close to the two-pion threshold has another ad- 
vantage: at these invariant-mass values the re-generation of sigmas from the pions will be small, 
i.e the generation of a cr-particle with this constrained invariant-mass, from pions produced by 
primary sigma decay, is statistically rare, which provides a justification for non-consideration of 
re-generation effects. Finally, the mean multiplicity of ex's, before reconstruction, within this win- 
dow is ~ 60. Ai*2(M) stands for the corrected factorial moments, which arise after the subtraction 
of the moments of the mixed events. The reconstruction procedure in the event-by-event analy- 
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FIG. 7: (color online)/?! graph a) we present the histogram of the exponents £ of second factorial moment, 
for all the event population, arising from a conventional moment analysis, in 5 — 5 case, for tf rz = 10 fm. 
In graph b) we show the corresponding values of R which determine the precision of the fits. In graph c) 
we depict the second factorial moment of a single event ( with slope less than one third standard deviations 
away from the histogram average), for tf rz = 10 fm. In graph d) we present the average second factorial 
moment, for tf rz = 10 fm. 



sis, compared to the reconstruction in the whole ensemble which is widely used in the literature, 
possesses some additional problems. Firstly, in a small fraction of events the subtraction of the 
moments of the mixed pions leads to negative moments. This is a sign that for these events, and 
at the specific cell partition, the reconstruction procedure did not work and therefore they have to 
be excluded. Secondly, another fraction of reconstructed events possesses positive moment values 
with slopes (i.e. (^-values) smaller that or greater than 1. This is also a result of an unsuccessful 
reconstruction and means that their profile is completely overwhelmed by the background of mixed 
events. Both these problem sources have their origin in the low statistics of the event-by-event 
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FIG. 8: (color online)/?! graph a) we present the histogram p(Q of the exponents £ of second factorial 
moment, for all the event population at the freeze-out (tf rz = 10 fm), performed after the reconstruction 
procedure with an analysis window of 10 MeV in the di-pion invariant mass, in 8 = 5 case. In graph b) we 
show the corresponding values of R 2 which quantify the linearity accuracy. In both a) and b) subfigures we 
have kept only those reconstructed events which possess R 2 > 0.8. In graph c) we depict the second factorial 
moment of a single reconstructed event (with slope less than one fifth standard deviations away from the 
histogram average), for the same universality class, analysis window and freeze-out time. In graph d) we 
present the average second factorial moment for the ensemble of reconstructed events. 



analysis, contrary to the ensemble averaged one. However, for the majority of events the procedure 
is successful, and this allows us to acquire the information which would be hidden in an average 
analysis. 

The time- evolution affects the initial geometry of the critical system in a two-fold way: Firstly, 
the power law weakens, i.e. most of the exponents £ are smaller than the critical value jfe = 2/3. 
Secondly, the linearity demanded in relation (16) is spoiled compared to an ideal fractal case. 
Both these features are obvious in fig. 8. Note that in order to ensure the power-law behavior we 
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cut in the final analysis the reconstructed events with Ft? < 0.8. These events either correspond 
to a partially unsuccessful reconstruction, or (which is the usual case) their time evolution has 
significantly spoiled the initial critical profile. In both cases they are of no interest for the purposes 
of this work, which focus on the survival of criticality. After their elimination we remain with 
about three fifths of the initial population. 

Although time evolution weakens the power-law strength and spoils its form, one can clearly 
observe strong signs of the initial criticality. That is, the slope £ ~ 0.42 of fig. 8d is rather large 
and linearity is satisfactory. This result is displayed more transparently in subfigure 8a. The 
peak of at small £ reveals the complete loss of the initial criticality in the majority of events. 
However, the structure of p(£) and especially its secondary peak is a robust sign that the initial 
fractal geometry has been survived in a subset of events. Quantitatively, we can calculate the 
average £ for those events of fig. 8a that are on the right of the minimum that separates the 
conventional background from the secondary peak. We find (C)|>min. ~ 0.58, which is very close to 
the initial value 2/3 of the critical index. Thus, evolutionary intermittency, in an event-by-event 
basis, provides a clear and unambiguous measure of revealing the signs of the survived criticality. 
Finally, we have confirmed that, although always carrying the information of the initial critical 
geometry, (C)|>min. is sensitive to the value of the isothermal critical exponent 5 (which determines 
the universality class) . This behavior enhances the robustness of evolutionary intermittency in the 
investigation of the QCD critical point, and furthermore it supports its use in order to trace the 
QCD universality class. 

V. SUMMARY AND CONCLUSIONS 

In this work we have explored, both in physical and momentum space, the dynamics of the 
critical isoscalar condensate which is expected to be formed near the QCD critical point. For this 
purpose we have adapted the <7-model Lagrangian in order to describe correctly the characteristics 
of the order parameter associated with the critical endpoint of the QCD phase transition [15] and 
we have developed a scheme to propagate the corresponding a-configurations in momentum space. 
The issue is of primary importance in the search for the existence and location of the QCD critical 
point, in experiments with nuclei. At the phenomenological level these fluctuations are expressed 
through the fractal mass dimension of the cr-field configurations, determining the properties of the 
condensate at criticality [7, 11]. 

In particular, we have considered the geometrical deformation of the transverse momentum 
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pattern of the critical cr's as a function of time, in terms of factorial moment analysis. To take 
into account dynamical effects related with this deformation we have introduced the concept of 
evolutionary intermittency, quantified by the total distribution of the critical indices £. We have 
calculated them by investigating the transverse momentum factorial moments in an event-by-event 
basis, assuming that each event is characterized by a different decay time lying in the freeze- 
out zone. Thus, evolutionary intermittency, describing event-by-event dipion pair correlations in 
transverse momentum space, is the basic observable in our approach. It incorporates, in terms of 
fractal mass dimension values, all the out-of-equilibrium dynamical history of the critical system. 
In this context, the main result of our treatment is the survival of critical fluctuations in momentum 
space even if dynamical effects prevail in the freeze-out era. This survival is exrpessed as a secondary 
peak in the distribution of the critical indices representing the evolutionary intermittency profile 
of the considered process. 

The description of the critical fluctuations for the isoscalar condensate within the presented 
approach is compatible, at least in an approximate scheme, with (a) the sigma dynamics, (b) the 
initial conditions of a second-order phase transition, (c) the effects in medium and in particular 
the partial restoration of chiral symmetry and (d) the requirements of observability in nuclear 
collisions. However, we do not pretend to present a complete picture since the introduced quench 
scenario does not allow for the parametrization of the freeze-out state of the formed critical 
condensate in terms of chemical potential and temperature, which are measurable quantities 
in a heavy- ion collision experiment. In addition we have neglected collective transverse flow 
phenomena. To incorporate these issues, one should include additional constraints in the dynamics 
which are beyond the present study. Last but not least, the assumption of pions originating 
exclusively from critical sigmas and not suffering final-state interactions, constitutes an idealized 
scenario for the evolution of the critical system. Nevertheless, despite its limitations, the treatment 
of the present work is a realistic step towards the determination of suitable observables associated 
with the critical fluctuations in a heavy- ion collision experiment. The obtained result are of 
interest especially in front of the forthcoming experiments at RHIC and SPS which are expected 
to reach this region in the QCD phase diagram. 
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